// STL
#include <cmath>               // std::sqrt()
#include <stdexcept>           // std::runtime_error()
#include <tuple>               // std::tie()

// jScience
#include "jScience/linalg.hpp" // Vector<>, Matrix<>, inverse(), dot_product()


//
// DESC: The Mahalanobis distance is a measure (of the distance) between EITHER
//
//     a point and a distribution           OR
//     two points of the same distribution.
//
// -----
//
// INPUT:
//          Vector<double> x   :: (first) point
//
//          Vector<double> mu  :: distribution mean (or second point)
//          Matrix<double> cov :: distribution covariance matrix
//
// OUTPUT:
//          double         d   :: distance measure
//
// =====
//
// NOTE: See:
//
//     https://en.wikipedia.org/wiki/Mahalanobis_distance
//
inline double distribution::Mahalanobis(
                                        const Vector<double> &x,
                                        // -----
                                        const Vector<double> &mu,
                                        const Matrix<double> &cov
                                        )
{
    double d;

    Vector<double> x_mu = x - mu;

    int            INFO0;   // DGETRF exit code
    int            INFO;    // DGETRI exit code
    Matrix<double> cov_inv; // Matrix<> inverse
    std::tie(
             INFO0,
             INFO,
             cov_inv
             ) = inverse(
                         cov
                         );

    if( (INFO0 != 0) || (INFO != 0) )
    {
        throw std::runtime_error( "error in distribution::Mahalanobis(): ( (INFO0 != 0) || (INFO != 0) ) on return from inverse()" );
    }

    d = std::sqrt( dot_product(x_mu, (cov_inv*x_mu)) );

    return d;
}
